%%% Implemenation: Anton Kidess <a.kidess@tudelft.nl>, TU Delft 2011

%%% Analytical solution of Shercliff's case
% Solution from p 47ff - Müller, U., Bühler, L., 2001. Magnetofluiddynamics in channels and containers. Springer. 
% URL http://www.worldcat.org/isbn/9783540412533
% and
% Shercliff, J. A., 1953. Steady motion of conducting fluids in pipes under transverse magnetic fields. Mathematical Proceedings of the Cambridge Philosophical Society 49 (01), 136-144. 
% URL http://dx.doi.org/10.1017/S0305004100028139

disp('Start');

%% settings
d = 1; %duct diameter
Ha = 20; % hartmann number

N = 50; %iterations on the fourier series
partitions = 501; %spatial discretization

write_vel = 1; %write solution to file
write_option = 'sets'; %write_option [surf|sets]
plot_solution = 0; %plot solution in matlab

%% initialization
lambda = zeros(N,1); %initialize lambda array (see ref.)
k = zeros(N,1); 
p1 = zeros(N,1); %pressure drop
p2 = zeros(N,1); %pressure drop
a1 = zeros(N,1);
a2 = zeros(N,1);

%set up coordinate mesh
z = linspace(-d,d,partitions);
y = linspace(-1,1,partitions);

%% Computation
disp('Setting up lambda, k, pi, ai');
for i = 1:N
    j=1+2*(i-1); %j is what the Mueller calls i
    lambda(i) = j * pi / (2*d);

    k(i) = 2 * sin(lambda(i)*d) / (lambda(i) * d);
    p1(i) = 0.5 * (Ha + sqrt(Ha*Ha + 4 * lambda(i) * lambda(i)));
    p2(i) = 0.5 * (Ha - sqrt(Ha*Ha + 4 * lambda(i) * lambda(i)));
    a1(i) = sinh(p1(i));
    a2(i) = sinh(p2(i));
end

disp('Computing fi, ui');
fiY = zeros(N,1);
f = zeros(N, partitions);
ui = zeros(N, partitions);
for i = 1:N
    fiY(i) = a2(i) * cosh(p1(i)) - a1(i) * cosh(p2(i));
    for j = 1:partitions
        f(i,j) = a2(i) * cosh(p1(i) * y(j)) - a1(i) * cosh(p2(i) * y(j));
        ui(i,j) = k(i) / (lambda(i) * lambda(i)) * (1 - f(i,j)/fiY(i));
    end
    %u(i,partitions) = 0;
end

disp('Reconstructing velocity field');
vel = zeros(partitions, partitions);
vmean = 0; %calculate mean flow
for j = 1:partitions %loop over y
    for k = 1:partitions %loop over z
        for i=1:N %loop over series
            vel(j,k) = vel(j,k) + ui(i,j) * cos(lambda(i) * z(k));
        end
        vmean = vmean + vel(j,k);
    end
end
vmean = vmean / (partitions * partitions)

%% plotting
if (plot_solution == 1)
    disp('Plotting solution');
    surf(z,y,vel/vmean);
    %contour(z,y,vel);
    xlabel('z');
    ylabel('y');
    axis('square');
end

%% Write to file
if (write_vel == 1)
    disp('Writing result to file');
    %write to file
    fid = fopen('shercliff_hr.txt', 'w');
        %header
        fprintf(fid, '#Shercliffs case\n');
        fprintf(fid, '#Ha = %d\n', Ha);
        fprintf(fid, '#vmean = %f\n', vmean);
        
        %% write surface
        if (strcmp(write_option,'surf'))
            disp('Write surface');
            fprintf(fid, '#z\ty\tvelocity\n');
            %data
            for j = 1:partitions %loop over y
                for k = 1:partitions %loop over z
                    fprintf(fid, '%6.5f %6.5f %6.5f\n', z(k), y(j), vel(j,k)/vmean);
                end
            end
        end
        
        %% write sets
        if (strcmp(write_option,'sets'))
            disp('Write sets');
            fprintf(fid, '#z\ty\tvelocity\n');
            j = round(partitions/2); %loop over y
            for k = 1:partitions %loop over z
                fprintf(fid, '%6.5f %6.5f %6.5f\n', z(k), y(j), vel(j,k)/vmean);
            end
            fprintf(fid, 'e\n\n\n');
            
            fprintf(fid, '#z\ty\tvelocity\n');
            k = round(partitions/2); %loop over y
            for j = 1:partitions %loop over z
                fprintf(fid, '%6.5f %6.5f %6.5f\n', z(k), y(j), vel(j,k)/vmean);
            end
        end
    %close file
    fclose(fid);
end
